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ABSTRACT 

High-resolution 2D-MHD numerical simulations have been carried out to in- 
vestigate the effects of continuing infall of clumpy gas in extended Hi galactic 
disks. Given a certain accretion rate, the response of the disk depends on its 
surface gas density and temperature. For Galactic conditions at a galactocen- 
tric distance of ~ 20 kpc, and for mass accretion rates consistent with current 
empirical and theoretical determinations in the Milky Way, the rain of compact 
high velocity clouds onto the disk can maintain transonic turbulent motions in 
the warm phase (~ 2500 K) of Hi. Hence, the Hi line width is expected to be 
~ 6.5 km s^^ for a gas layer at 2500 K, if infall were the only mechanism of 
driving turbulence. Some statistical properties of the resulting forcing flow are 
shown in this Letter. The radial dependence of the gas velocity dispersion is also 
discussed. 

Subject headings: galaxies: intergalactic medium — galaxies: ISM — hydrody- 
namics — ISM: kinematics and dynamics — ISM: structure — turbulence 



1. Introduction 



Hi line widths, aw, are observed to vary from ~ 12 to 15 km s~^ in the central parts 
to 6-8 km s~^ in the outer parts. Beyond the optical disk, the maps of aw display a 
patchy distribution with values populating the interval from 5 to 9 km s~^, independent of 
galactocentric distance. The median value is fairly constant along the extended outer parts 



with about the same univ ersal value 7 ± 1 km s ^ for all the observed galaxies (iLewis 



1984 



Dib. Bell fc BurkertI l2006l . and references therein). By fitting two-component Gaussians to 
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the Hi profiles in NGC 6822, the second-moments for the cold an d warm phases of the H I 
were found to be ~ 5 km s~^ and oy^ ~ 8 km s~^, respectively (Ide Blok fc Walterll2006l ). 



The observed are likely to represent turbulent small-scale motions: the existence of a 
cold phase in the outer disks, the observed level of star formation and the double exponential 
radial profile of the star formation rate have been also interpreted as a consequence of the tur 
bulence gas compressions in the outer disks (e.g., 
2006h . 



Ferguson et al.lll998l ; lElmegreen fc Hunter 



Dib. Bell &: BurkertI (120061 ) found that if turbulence is driven by supernova, the velocity 
dispersion of the Hi gas in the quiescent regime is ~ 3 km s~^, a factor of ~ 2 smaller 
than the observed values. Thus, either there are other physical processes driving turbu- 
lence or the supernova feedback efficiency has been underestimated. There exist several 
physical mechanisms that could drive turbulence even in the absence of star formation: hy- 
drodynamic or magnetohydrodynamic instabilities, frequent minor mergers of small satellite 
clumps, ram pressure or infalling gas clouds. Amongst t he hydro dynamical instabilities, the 



thermal instability by itself cannot sust ain turbulence (jSanchez-Salcedd 1200 ll : iGazol et al. 



200ll : iBrandenburg. Korpi fc Med 120071 ). The magnetorotational instability, on the other 
hand, can only account for an amount of turbulent motions ~ 4.5 km s~^ for the three 
components, if thermal broadening is not subtracted, and might be completely suppressed 
by stellar feedback (Dib et al. 2006 and references therein). 

Our main goal is to assess how much of the velocity dispersion observed in the ISM 
can be due to the impact of high velocity clouds (HVCs) and intermediate velocity clouds 
(IVCs). Ample evide nce for the pres e nce o f continuing gaseous infall to the Galactic disk 
has been compiled in iBeckman et al.l (120031 ). All the studies firmly suggest that the infiow 
in the Galactic disk has been constant, with an accretion rate of ~ 1 M© yr~^, or has even 
increased, during its lifetime. H I observations have revealed accretion of both diffuse and 



discrete structures from the extended en vironment in M31 and M33 (IThilker et al.l 12004 



Braun &: Thilkerl |2004| ). iPutmanl (120061 ) infers that the present total mass in condensed 
infalling clouds around the Milky Way is ~ 6 x 10® Mq if they are all at distances < 60 
kpc. HVCs may be a repository for large amounts of gas if clouds fall into the disk rap idly 
after they are formed jMalonev fc Putman|[2003l : iMaller fc BuUockl 120041 : |Putman|[2006h . In 
the next sections, we consider the turbulent Hi dynamics driven by the infall of clouds. 
Preliminary results were reported in Sanchez- Salcedo et al. (2007). 



- 3 - 



2. The model 

The ideal MHD equations are integrated using the ZEUS code (Stone & Norman 
1992a,b). A local Cartesian frame of reference with x and z corresponding to the hori- 
zontal and vertical directions, respectively, was adopted to simulate a small patch of our 
Galaxy at a distance Rq from the Galactic center. The galactic symmetry plane, z = 0, is 
placed in the middle of the computational domain. The simulation domain is a square of 
size L X L with 1024^ zones. The ambient medium that is interacting with the HVCs is 
initially in plane-parallel magnetohydrostatic equilibrium with scaleheight ho in an external 
gravitational potential. Here the scaleheight is defined as half width at half maximum of 
the vertical volume density of the gas layer. The ambient gas is initially isothermal in space 
with thermal sound speed Cs and its evolution is nearly isothermal with specific heat ratio 
7 = 1.01. Gas self-gravity was disregarded. The magnetic field has only one component 
along the |/-axis and is initially stratified in the ^-direction, By{z). Thus, there is no mag- 
netic tension, only magnetic compression. Two magnetic configurations are explored. The 
first one is the high-latitude e xtension of the Galact ic thick disk of Boulares & Cox (1990, 



BC model hereafter) following ISantillan et al.l (119991 ). The vertical gravitational force and 



the magnetic field configuration were scaled to the outer Galaxy, neglecting edge effects. 
The second set of models assumes that By falls off with the gas density (Spitzer case) in 
the initial equilibrium configuration, so that a, defined as the ratio between the magnetic 
and thermal pressure, is constant. In all the simulations, the size of the domain is L > IO/iq 
in order to avoid spurious boundary effects. We apply periodic boundary conditions in the 
X-axis, and open boundary conditions on the upper and lower faces of the computational 
domain. 

The too low thermal gas pressure at the outer parts of ext ended galactic disks will permi t 



only the warm (several times 10^ K) phase of H I to exist (e.g.. lElmegreen fc Parravandll994j ). 



Heiled ( 120011 ) provides observational evidence that about half of the mass of the diffuse gas 
in the Galaxy may have temperatures from a few hundred to a few thousand kelvin. In this 
work, we will explore values in the range 4.5-8 km s~^. 

The ambient gas is subject to a continuous rain of HVCs. In most simulations, if not 
specified otherwise, all the clouds with initial radius Rd and uniform internal density nd, 
are injected at z = Zd with a vertical velocity Vd, and are placed randomly in x. We found 
that the statistical properties of the flow in the disk are similar when the injection occurs 
along both caps z = ±Zd- 

The remaining five dimensionless parameters that characterize our simulation models 
are a{z), Rd/ho, Vd/cg, Ud/no, and Aacc, where uq is the midplane volume density at t = 
and Aacc is the fraction of accreted mass in one sound crossing time, Aacc = {ho/cg) (S/Sj) 
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with Sj is the initial column density of the disk. Note that once Cg is fixed, the scaleheight ho 
is determined by the depth of the external gravitational potential. The collisions of infalling 
clouds onto the gaseous disk will take the gas out of hydrostatic equilibrium. Consequently, 
the scaleheight h of the H I layer may evolve in time. 

3. Results 

The disk gas is subject to a stocastic driving force caused by the MHD pull of randomly- 
injected clouds, plus the restoring gravitational force which tends to push back the gas to 
the midplane. The level of substructure in the resulting fiuid depends on the structural 
parameters of the clouds and injection velocity. We focus first on describing the evolution 
of our fiducial BC model, labelled as run BCa, that is intended to represent conditions at a 
galactocentric distance Rq = 20 kpc. In this model, the parameters of the disk are Cg = 8 km 
s~^, E ~ 1.2 Mq pc~^ and ho = 220 pc. For the clouds we adopted Vd = 100 km s~^, Zd = 2 
kpc, rici = 0.1 cm~"^ and Rd = 50 pc (see Sanchez-Salcedo et al. 2007 for a discussion of this 
choice). The injection rate of clouds in the whole box domain is 300 per Gyr, corresponding 
to an accretion rate defined as AttRIJ] of 0.6 M0 yr~^. While our simulations are 2D, this 
situation is energetically equivalent to an accretion rate of 750 clouds per Gyr and per kpc^ 
of galactic disk. Since the kinetic energy of a cloud is ~ 10^° erg, the energy input rate 
from impacting clouds is 7.5 x 10^^ erg Gyr~^ kpc~^, which is equivalent to 75e~^ supernova 
explosions per Gyr and kpc^ of disk, each with mechanical energy e x 10^^ erg, where e ~ 0.1 
is the efficiency factor. At the end of the simulation (2 Gyr), the total mass in our grid 
increased by a factor 1.45. 

Clouds interact with the ambient medium and dilute (see FigH]). At a height of 500 pc 
they lose their identity and mix up with the turbulent disk. Hence the number of 'active' 
clouds is 300 Gyr^^ x (1.5 kpc/100 km s~^) ~ 4-5. The infalling clumpy fiow perturbs the 
disk through MHD waves that interact in a complicated fashion, producing a network of 
plumes and shells (see Fig. [T]). As a consequence, the ID rms velocity of the gas, a, increases 
in time and saturates after 100 Myr, reaching values that oscillate between 6 and 7.5 km s~^ 
around the mean value a (see Fig. [2]). The difference between ax and is only ~ 0.5 km s"^. 
The velocity dispersion in each direction was obtained by fitting the mass-weighted velocity 
profile for cells within \z\ < hf, with hf the asymptotic value of the scaleheight. Examples 
of mass-weighted velocity profiles are given in Figure [2)d,c. Remind that aw = (c"^ + c^Y^"^ 
for isothermal gas, thus, aw ~ 10 km s~^. For Vd = 50 km s~^, ax and az reach similar 
values as in the case Vd = 100 km s"^. These results suggest that continuous accretion of a 
clumpy gas may contribute significantly to the random motions observed. 
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The vertical density profile averaged on horizontal cuts is shown in Fig.[3l Note that even 
though the injection of clouds occurs only through the upper cap z = z^i, the density profile 
is quite symmetric and the net vertical displacement is relatively small. The scaleheight is 
~ 640 pc at ~ 2 Gyr, similar to the observed value at i? = 20 kpc. 

The injected energy is mostly dissipated in radiative shell compressions and behind shock 
fronts, being the dissipation timescale for turbulent kinetic energy ~ 8 Myr. The energy 
spectrum contains information about how the kinetic energy is distributed at different scales 
and about the nature of the driving mechanisms of the random motions in the compressible 
ISM. Due to the stratification and the anisotropy of the random motions, we have calculated 
the power spectrum along horizontal cuts for each component of the velocity separately. The 
power spectrum follows a power-law with different exponent for each component Ei{k) ~ 
/c~°\ with ttr = 3 and a. = 4. 



Vazquez-Semadenil (Il994j ) and iPadoan. Nordlund fc Joned (119971 ) reported lognormal 



density probability distribution functions (PDF) for two- and three-dimensional isothermal 
turbulent flows. We have checked that, in fact, our PDF is approximately lognormal with 
excesses at large and low densities probably because in our case the turbulence is not isotropic 
and, in addition, there is some probability of finding density substructure that has not a 
turbulent origin, it being the remnants of infalling clumps. 

When the clouds are randomly injected along both caps, the mass enhancement in the 
central disk {\z\ < 600 pc) is slightly larger (about 10%) than when the injection only 
proceeds along one side, because the symmetry in the injection leads to a more efficient 
balance of momentum. However, a{t) is quantitatively similar and thus there is no need to 
show it again. 

The elasticity of the disk to collisions depends on the speed of sound and on the magnetic 
field strength. In order to see the sensitivity of a on the adopted magnetic configuration, we 
started with an a = 1 Spitzer model and the same Cg and total pressure at the midplane as 
in run BCa. ax and (Xz were found to be 5 and 6.2 km s~^, respectively. 

It turns out very difficult to infer from basic principles how a scales with Cg or with 
any other parameter (Sanchez-Salcedo et al. 2007). In order to gain some insight, Spitzer 
configurations with a = 1 but different Cg, between 5 and 12 km s~^ (spaced by 1 km s~^), 
were considered. If the total pressure at the midplane 2; = is identical in all the models 
and equal to the value in the standard modeQ, we find that a scales linearly with c^; a good 



The fact that the total pressure at z = is the same in these models implies that the surface gas density 
increases for low values of Cs ■ 
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fit to the data is cTz = 0.55cs + 1.16 (thermal broadening has been subtracted) when Vd = 50 
km s~^. If, instead of the total pressure, the surface gas density is taken the same as in the 
case BCa in all the experiments, then oc y/c^. For instance, using again Vd = 50 km s~^, 
(j^ varies from 7.4 km s~^ to 5.0 km s~^ when the sound speed drops by a factor of 2 (from 
12 km s^^ to 6 km s~^). 

In order to simulate various galactic conditions at different galactocentric distances, 
runs have been carried out varying Aacc and Rd/ho. Assuming that Rd does not depend 
strongly on radius, Rd/ho should decrease with R because the observed h increases with R, 
but it is uncertain as Rd could also increase with R as well. The radial variation of Aacc is 
also uncertain because the accreting mass flux E(i?) and the temperature structure Cs{R) 
should be known. As a departure assumption consider that accretion of an intergalactic 
plane-parallel flow at inflnity occurs ballistically S oc R~^/'^ (Lopez-Corredoira et al. 2002). 
For our Galaxy, S decreases by a factor of ~ 10 between 15 and 25 kpc, whereas h increases 
by a factor of 4. Thus, Aacc varies by a factor of ~ 30 assuming that Cg is constant with R. 
In order to bracket conditions at these galactocentric distances, Aacc has been varied by a 
factor of 25 for run BCb as compared to run BCc, whereas Rd/ho varies by a factor of 4 (see 
Tabled]). In runs BCb and Sp-b, reaches a maximum and decreases afterwards partly 
because the accreted mass in the disk is no longer negligible in the length of the run. Since in 
these experiments Vz deviates signiflcantly from a Gaussian distribution, whereas Vx is still 
well described by a Gaussian, only the values of a^/cg are reported in Tabled! According to 
run BCc, which represents Galactic conditions at i? ~ 15 kpc, we flnd that for Cg = 4.5 km 
s~^, ^ 3.4 km s~^ and aw = 5.6 km s~^ if thermal broadening is taken into account. At 
Cs = 8 km s~^, ax ~ 4.5 km s""*^ (thus aw = 9 km s""*^). At the outer far disk (~ 25 kpc), ax 
is found to be 7.5-10 km s~^ for Cs = 4.5 and 8 km s~\ respectively {aw ~ 9-13 km s^^). 
To derive these values, we rescaled the results of run BCb appropriately using our previous 
result that ax oc y/c^. Slightly smaller ex's are derived for Spitzer models (see Tabled]). Line 
widths of ~ 10 km s~^ as predicted at the edge of the disk seem rather high, but note that 
empirical deter nainations are scarce and very u ncertain at low surface gas densities (S < 1 



M© pc-2) fe.g.- IShostak fc van der Kruitlll984] ). 



If accretion occurs through rare mergers of very massive clouds, it would be hard to 
see how such a localized events could produce the observed uniform level of turbulence. 
Nevertheless, experiments with input rates of new clouds 25 times smaller than in our flducial 
simulation (i.e., 12 per Gyr in our box domain) but cloud radii 5 times larger (250 pc in 
physical units), so that the mass accretion rate is unchanged, still maintain a rather uniform 
level of turbulence everywhere. 
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Discussion and conclusions 



We have shown that rms velocities of 3-8 km s^^ in the warm phase of Hi naturally 
arise if outer disks are continuously stirred by a clumpy accretion flow at a rate consistent 
with observations. The exact radial variation of the H I line width depends on the uncertain 
contribution of the thermal broadening as a function of radius and on the adopted accreting 
mass flux which ultimately depends on which the origin of the condensing clouds 

is. An interesting possibility is that the outer disk is being ra ined on by Galactic fountain 
material driven by inner starbursts (IBenjamin fc Shapirdll993l ) or due to matter circulation 
in the halo. Under the rather extreme assumption that clouds fall ballistically from infinity, 
a\Y is expected to vary only ~ 60% between 1 and 1.7 times the optical radius. The coupling 
between gas accretion and stellar feedback could account for the uniformity of the Hi line 
widths. 

So far, our quasi-isothermal simulations do not capture all the physics of the multi- 
phase ISM to permit a detailed comparison between predictions and observations. Three- 
dimensional simulations including cooling and heating in thermally bistable medium, plus 
a spectrum of mass and radius of the clouds are being undertaken. These experiments will 
allow us to derive the velocity dispersion for each component (cool, warm and hot phases) 
and compare them with high-resolution Hi observations that have bee n able to identify 
cool and warm neutral components in the ISM of external galaxies (e.g., Ide Blok fc Walter 
20061 ). The accretion of clumpy gas can induce turbulence-compressed regions and trigger 
star formation, which should be also included. 



It is important to consider other phenomena arising in the present accretion scenario. 
For instance, if turbulence is driven by oblique accretion flows, a large asymmetry in the 
scaleheight of the H I disk as that observ ed between the northern and southern halves of the 
Milky Way (ILevine. Blitz fc Heilesll2006l) is expected because oblique flows have an azimuthal 
dependence (see ISanchez-Salcedd 120061 ). In addition, infalling clouds moving supersonically 
at the upper warm galactic disk develop optical-emitting shocked regions. At a density of 
the preshocked gas of 5 x 10~^ cm~^, the face-on Ha surface brightness {Ia)± of a certain 
shocked region m ay range betwee n 10 and 32.5 mR for shock velocities of 50 and 100 km 



s ^, respectively (IRaymondl Il979l ) . In an edge-on galaxy, the expected surface brightness. 



on scales of ~ 1.5 kpc, is the result of the contribution of all the shocked regions through 
the sight line: I ^ {Sn /2) NdRliR^Ia) ±, where N^i is the number of clouds per unit volume, 
in our models Nd ~ 2.5-5 kpc~^. At -R = 20 kpc, I = 12-40 mR depending on the shock 
velocity, implying emission measures EM = J n^dl = 0.03-0.1 pc cm~^, and assuming gas 
at 8000 K. While the detection of such a Ha brightness is challenging, extremely deep Ha 
images of edge-on external galaxies might constrain the form in which galaxies accrete mass. 
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Table 1. Some relevant runs 



Run" 


Vd/Cs 


rid/no 


Rd/ho 


10' A,^, 


C^x/Cs 


BCa 


12.5 


1 


0.23 


1.2 


0.75 


BCb 


12.5 


1 


0.11 


6.0 


1.25 


BCc 


12.5 


1 


0.46 


0.24 


0.56 


BCd 


12.5 


0.1 


0.23 


1.2 


0.25 


BCc 


1.25 


1 


0.23 


1.2 


0.69 


Sp-a 


8.33 


0.9 


0.17 


2.2 


0.70 


Sp-b 


8.33 


0.9 


0.17 


11.55 


1.05 



^ BC refers to the extension of the Boucoulares & Cox (1990) 
model. Sp indicates that a Spitzer model with a = 1 was used. 
' In physical units, an accretion rate of 0.6 M© yr~^ corresponds 
to 1.2 in the BC model and to 2.2 in the Spitzer model, for the 
Galactic surface density at 20 kpc. 
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X (kpc) 

Fig. 1. — Density {color logarithmic scale) and velocity field (arrows) at t = 1.7 Gyr for our 
run BCa. A general view of the disk between —2 kpc < z < 2 kpc (note that the 2;-axis has 
been cut-off) is shown (top). A magnification zoom right into the marked region is displayed 
[bottom). Notice the presence of vortical motions and the dilution of the infalling clouds. 
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Fig. 2. — Horizontal {solid line) and vertical {dashed line) velocity dispersions, with thermal 
broadening subtracted, for our fiducial case (run BCa), as a function of time (panel a). 
Panels b and c show the mass- weighted profile at i = 1.9 Gyr. The profiles sample those gas 
elements within one scaleheight of the disk. Overplotted are the Gaussian fits (crosses) . 
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Fig. 3. — The vertical density profile averaged on horizontal cuts is shown at three selected 
times, t ~ 0, 1.27 and 1.9 Gyr, for run BCa. The corresponding scale heights are 220, 660 
and 640 pc, respectively. 



